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ABSTRACT 


We  consider  the  effects  of  small 
random  perturbations  on  deterministic 
systems  of  differential  equations.  The 
deterministic  systems  of  interest  have 
oscillatory  dynamics  and  may  undergo  a 
bifurcation  (the  Hopf  bifurc^ti.9^.  -We-4 
^formulate  a first  exit  probiem^f oV  ex- 
periments beginning  near  stable  and 
unstable  limit  cycles.  The  unstable 
limit  cycle  is  surrounded  by  an  annulus. 

Of  interest  is  the  probability  of  first 
exit  from  the  annulus  through  a specified 
boundary,  conditioned  on  initial  position. 
The  diffusion  approximation  is  used,  so 
that  the  conditional  probability  satisfies 
a backward  diffusion  equation.  Appropriate 
solutions  on  the  backward  equation  are 
constructed  by  an  asymptotic  method.  The 
behavior  of  the  stochastic  system  in  the 
vicinity  of  stable  and  unstable  limit 
cycles  is  compared.  When  the  deterministic 
system  exhitits  the  Hopf  bifurcation,  the 
above  analysis  musf"Ee  modified.^. 

Uniform  solutions  of  the  backwarcT 
equation  are  constructed.  The  solutions 
are  analogous  to  Hadamard's  solution  of 
the  point  source  problem  for  the  wave 
equation.  Numerical  examples  are  used 
to  compare  the  theory  with  Monte  Carlo 
experiments. 
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INTRODUCTION 


In  the  past  few  years,  the  anlaysis  of  oscillatory  non- 
linear dynamical  systems  has  received  considerable  attention, 
due  to  a variety  of  biological,  physical,  and  chemical  applications. 
Mainly,  the  analysis  has  been  based  on  systems  of  deterministic 
differential  equations  and  involved  classification  of  dynamical 
behavior  of  the  systems.  Most  of  the  analyses  ignored  fluctuations 
that  are  always  present  in  such  systems  (Ludwig, (1975) , Van  Kampen 
(1976) , and  White  (1977)  are  exceptions) . In  this  paper,  we 
consider  the  effects  of  fluctuations  on  systems  with  oscillatory 
behavior.  We  consider  an  autonomous  system 


x = b(x)  xeR2 


(1.1) 


that  has  a periodic  solution (s).  Three  types  of  periodic  solutions 
are  of  interest  here:  l)a  fixed,  stable  limit  cycle,  surrounding 
an  unstable  focus  (figure  1A) ; 2)a  fixed  unstable  limit  cycle, 
surrounding  a stable  focus  and  enclosed  by  a fixed,  stable  limit 
cycle  (figure  IB) ; 3)the  Hopf  bifurcation  problem:  the  deterministic 
dynamics  depend  upon  a parameter  y.  As  ylO,  a stable  limt  cycle 
coalesces  with  an  unstable  focus  (at  y=0) . The  limit  cycle  dis- 
appears and  the  focus  becomes  stable  (figure  1C) . A "dual" 
bifurcation,  in  which  an  unstable  cycle  and  stable  focus  coalesce, 
is  shown  in  figure  ID  . 
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The  three  types  of  oscillatory  solutions  arise,  for  example, 
in  theoretical  population  dynamics  (e.g.,  Bazekin,  1975)  and  chemical 
reaction  dynamics  (e.g.,  Cohen,  1972).  The  fixed  unstable  limit 
cycle  occurs  in  the  treatment  of  molecule  ion-molecule  collisions. 

The  stable  limit  cycle,  with  superposed  fluctuations,  was 
treated  by  Ludwig,  White,  and  Van  Kampen.  We  include  it  here  for 
two  reasons.  First,  our  treatment  is  slightly  different  from  the 
others.  Second,  it  is  interesting  to  compare  the  stochastic 
dynamics  of  stable  and  unstable  limit  cycles.  The  unstable  cycle 
contained  by  a stable  cycle  arises  in  chemical  dynamics  (Tyson,  1977, 
Uppal  et  al,  1974).  In  the  engineering  literature,  the  unstable 
limit  cycle  is  called  a "hard"  oscillation,  and  the  stable  limit 
cycle  a "soft"  oscillation.  The  Hopf  bifurcation,  and  dual  Hopf 
bifurcation,  arise  in  many  situations  (Marsden  and  McCracken,  1976) . 
Our  interest  is  again  motivated  by  chemical  reaction  dynamics 
(Cohen,  1972,  Uppal  et  al,  1974).  Fluctuation  effects  have  not  been 
considered  in  these  systems. 

When  fluctuations  are  superposed  upon  the  deterministic 
dynamics  (1.1),  a number  of  interesting  questions  arise.  The 
type  of  question  that  should  be  posed  depends  upon  the  type  of 
deterministic  dynamics,  as  is  to  be  expected.  First,  consider  the 
stable  limit  cycle  (figure  1A).  Since  the  deterministic  attraction 
is  always  towards  L , the  question  of  interest  involves  how  fluc- 
tuations may  derive  the  system  away  from  L . Let  x(t)  denote 
the  random  variable  obtained  by  superposing  fluctuations  on  (1.1). 


In  this  case,  x(t)  in  (1.1)  is  the  appropriate  conditional  average 

% . 

of  x(t).  Let: 

v(t,  x)dx  = Pr{x  s x(t)  s x + dx}  (1.2) 

Thus  v ( t , x ) is  the  probability  density  for  x(t) . It  represents 
a natural  function  describing  the  stochastic  dynamical  system 
obtained  from  figure  1A.  If  we  let  t -*  °°  , then  v(t,x)  -»  v(x), 
the  equilibrium  or  stationary  density  of  eventually  finding  the 
process  between  {x,  x + dx}.  We  note  that  v(t,x)  is  independent 
of  x(0).  Namely,  we  are  solely  interested  in  the  forward  time 
density,  given  some  initial  distribution,  vQ(x(0)). 

On  the  other  hand,  the  initial  point  is  crucial  when  con- 
sidering an  unstable  limit  cycle  U (figure  2) . A phase  point 
initially  in  the  vicinity  of  U leaves  any  neighborhood  of  U 
with  probability  1.  Even  if  x(0)€  U,  fluctuations  will  drive  the 
point  away  from  U.  Ideally,  we  would  like  to  calculate  the 
probability  that,  given  x(0)  = x,  the  phase  point  reaches  a 
neighborhood  of  the  node  P rather  than  a neighborhood  of  the 
stable  limit  cycle  L.  In  general,  this  problem  is  too  difficult 
to  solve.  We  can,  however,  analyze  the  following  problem.  Let  s 
measure  distance  normal  to  U . Consider  two  contours: 


S 


1 


{x:  s(x)  = s1} 


{x:  s(x)  = s 2 } 


(1.3) 


(see  figure  2)  , with  s-^  > 0,  s2  < 0.  Consider  the  probability 

A/ 

u(x,t'  = Pr{by  time  t,  x(x)  has  left  the  annulus  (Sj_,S2), 
through  S2|x(0)  = x} 

This  probability  is  dependent  solely  upon  initial  position.  The 
stationary  equivalent  of  (1.4)  is  u(x),  which  is  the  probability 
that  x(t)  first  exits  from  (S^,S2)  through  S2> 

For  the  case  of  a system  exhibiting  the  Hopf  bifurcation, 

rw 

(figure  1-C) , we  are  again  interested  in  the  density  for  x(t). 

Now  the  density  v = v(x,t;y),  where  y is  the  parameter  charac- 
terizing the  deterministic  bifurcation.  Now  consider  the  dual  Hopf 
system  (figure  1-D) . For  small  y , a phase  point  will  leave  a 
neighborhood  of  P or  U and  approach  L with  probability  1. 

The  singularity  at  P/U  for  y=0  will  be  evidenced  by  very  slow 
deterministic  repulsion  from  P.  Let  L be  a neighborhood  of  L and 

T (x ) = E{  t : X(t)€  L , x(s)^ij,  s<t|x(0)  = x,  x(t)  reaches  L} 

Thus,  T (x)  is  the  expected  time  to  reach  L,  given  that  x(0)  = x. 

In  order  to  calculate  the  above  quantities,  we  need  to 
introduce  a stochastic  kinetic  equation.  In  section  2,  we  first 
specify  the  deterministic  dynamical,  equations  corresponding  to  the 
systems  pictured  in  figure  1. 


Next,  we  modify  these  equations  by  the  addition  of  a random  function. 
This  is  the  Langevin  approach.  We  obtain  a stochastic  kinetic 
equation  that  is,  usually,  too  difficult  to  treat  directly.  We 
treat  the  kinetic  equation  by  the  diffusion  approximation  of 
Papanicolaou  and  Kohler  (1974).  in  this  approximation,  v(x,t), 
u(x,t)  and  T(x)  all  satisfy  partial  differential  equations.  A 
small  parameter,  characterizing  the  intensity  of  the  fluctuations, 
arises  in  derivation  of  the  stochastic  equations. 

In  section  3,  we  analyze  canonical  problems  corresponding 
to  stable  and  unstable  limit  cycles  and  the  Hopf  bifurcations. 

Various  incomplete  special  functions  arise  in  the  analysis  of 
these  canonical  problems.  These  functions  are  generalized  in 
section  4,  where  we  calculate  v(t,x)  and  u(t,x)  by  the  use  of 
formal  asymptotic  methods,  for  stable  and  unstable  (fixed)  limit 
cycles.  The  stationary  solutions  v(x),  u(x)  have  interesting 
interpretations  in  terms  of  "hindsight"  and  foresight."  In  section 
5,  we  construct  v(t,  x)  and  u(t,x)  for  Hopf-type  dynamical 
systems.  We  show  that  the  solutions  in  section  4 breakdown  and 
how  the  uniformly  valid  solutions  can  be  obtained.  In  section  6, 
we  present  some  numerical  solutions  indicating  the  phenomena 
discussed  in  sections  4-6. 
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SECTION  2 


DETERMINISTIC  AND  STOCHASTIC  i.INETIC  EQUATIONS 


We  first  characterize  the  deterministic  equations  that 
lead  to  the  phase  portraits  of  interest.  Then  we  formulate  the 
stochastic  kinetic  equations  and  the  diffusion  approximation. 


DETERMINISTIC  KINETIC  EQUATIONS 
We  assume  that 

x^  = b:*'(x,y)  x(ER2  n€R^ 

has  a periodic  solution  $ . Let  P be  a steady  state  of  (2.1) 
b1(P,y)  = 0 for  all  i . Introduce  new  coordinates:  s,  which 
measures  distance  normal  to  4>  , and  0 which  measures  distance 
along  $ . (If  $ were  a circle,  then  s = r and  0=6). 

From  (2.1),  we  obtain  equations  for  s and  0 


s = b (s , 0 ,y  ) 
0 = b°  (s, 0 ,y ) 


The  variable  0 is  periodic,  with  period 
The  limit  cycle  is  stable  if 


(bs (0, 0 , y) ) < 0 
3s 


and  is  unstable  if 


|i(bs(O,0,y  ) ) > 0 


(2.1) 


(2.2) 


(2.3) 


(2.4) 


i 


If  9/9 s (bS(0,6,y))  - 0,  then  the  limit  cycle  is  neutrally  stable 
(or  unstable) . This  phenomenon  occurs  at  the  Hopf  bifurcation. 

To  consider  the  Hopf  bifurcation,  we  return  to  (2.1). 

Let  B = (bV  (p , y ) ) , and  let  X(y),  X*(y)  denote  the  eigenvalues 
of  B.  The  Hopf  bifurcation  is  characterized  by  the  following 


1)  When  y < 0,  X(y)/X*(y)  are  located  in  the  left-half 

plane . 

2)  When  y = 0,  the  eigenvalues  are  located  on  the  imaginary 
axis . 

3)  When  y > 0,  the  eigenvalues  are  located  in  the  right- 
half  plane.  Also,  the  condition 

d 


dy 


Re  X ( y ) 


= Yi  t 0 


y=0 


holds.  There  are  analogous  (dual)  conditions  for  the 

dual  Hopf  bifurcation.  Let 

i0  1 , . 2 

z = re  = x + ix 

Fenichel  (1975)  (see  also  Arnold,  1972)  has  shown  that  (2.1)  can 
be  put  into  the  form  (for  small  y) 


(2.5) 


(2.6) 


• 3 it 

r = ± (b^  - ny  r)  = b (r,c|>,ri) 


= X 2 + b2r2  + ny2r  = b<^(r , <t> , n) 


(2.7) 


i 


where  r = r(s,6),  <J>  = <J>(s,0)  are  regular  functions,  is 

defined  in  (2.5),  X2  > ° and  blfb2  ? 0 . The  (±)  sign  in 
(2.7)  is  included  so  that  both  the  Hopf  bifurcation  and  dual  Hopf 
bifurcation  can  be  treated.  The  function  n = n(y)  is  regular 
and 

n (0)  = o 


(2.8) 


At  the  bifurcation,  n = 0,  we  note  that 


br  (0 , $ , 0)  = bf^  (0, 4> , 0)  = bfrr  (0,$,0)  = 0 

(0 , 4> , 0)  / o . 


b^ 


rrr 


These  conditions  will  be  used  later  (sections  3 and  4). 

STOCHASTIC  KINETIC  EQUATION  AND  DIFFUSION  APPROXIMATION 

Equation  (2.1)  is  approximate  in  that  it  completely  ignores 
fluctuations.  On  the  other  hand,  deterministic  equations  (e.g., 
the  "law  of  mass  action"  in  chemistry)  often  yield  correct 

predictions.  Such  deterministic  equations  are  successful 
for  two  reasons.  First,  the  fluctuations  are  of  small  intensity. 
Second,  the  fluctuations  occur  on  a time  scale  rapid  compared  to 
the  macroscopic  equation.  Accordingly,  we  replace  (2.1)  by  a 
Langevin-like  equation  for  the  random  variable  xa(t): 


dx 

a 

~SE 


= b1(xQ(,y)  + f^  (xa)Yj  (t/a2)  . 


(2.9) 


(2.10) 
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In  (2.10),  (s)  is  a stationary,  zero  mean  process,  satifying 

the  mixing  condition  of  Papanicolaou  and  Kohler  (1974).  The 
parameter  e,  0 < e«l,  characterizes  the  intensity  of  the 
fluctuations.  In  chemical  systems 
E - Ve/V 

where  V is  the  volume  of  the  reacting  system  and  Vg  is  the 
elementary  volume,  i.e.,  the  volume  of  a sub-unit  of  the  reacting 
system  (Kubo  et  al,  1973,  Van  Kampen,  1976).  The  parameter  a, 

0 s a « 1 characterizes  the  time  scale  of  the  fluctuations.  As 
a -*  0,  xa(t)  converges  to  a diffusion  (Papanicolaou  and  Kohler, 
1974).  The  density  v(t,x)  defined  in  section  1,  satisfies 
(Papanicolaou  and  Kohler,  (1974)): 

vt  = I^al^v)ij  " * (t)i  + eci>v>i  . 

The  probability  u(t,x)  satisfies 

e ij  . /t_i 

ufc  = 2s  uij  + (b  + ec  )u^. 

The  expected  time,  T(x)  , satisfies 
-1  = + (b1  + ec1)Ti 

In  (2.12  - 2.14),  subscripts  indicate  partial  derivatives  and 
repeated  indices  are  summed  from  1 to  n.  Also, 

a1^  (x)  = fkfi(Ykl  + Ylk) 
ci(x)  = yklfj  -K  f* 

K 9X3  1 


(2.11) 


(2.12) 


(2.13) 


(2.14) 


(2.15) 
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where 

Ykl  = f E[Yk  (s)Y1  (0)]  ds 
0 


In  later  sections,  we  will  specify  the  appropriate  boundary  condi- 
tions for  (2.12  - 2.13).  Equation  (2.14)  is  treated  elsewhere 
(Mangel , 1977 ) . 


SECTION  3 


CANONICAL  PROBELMS  AND  SPECIAL  FUNCTIONS 


We  now  consider  xCR1  , c1  = 0.  The  stationary  solutions 
versions  of  (2.12,  2.13)  and  appropriate  boundary  conditions  are 

00 

0 = f(av)xx  “ (bv)x  ; f v(x)dx  = 1 ;▼-*■<>  as  |x|-*-~  (3.1) 


0 = |au  + bu  u(x.)  = 0 u(x~)  =1  x.  < x_  (3.2) 

2 xx  x 1 2 12 


We  assume  that  the  deterministic  system 


x = b (x,  y) 


(3.3) 


has  a steady  state  at  x = xQ.  The  deterministic  equation  (3.3) 
corresponds  to  the  following  degenerate  "planar"  dynamical  system. 
Let  (r,9)  denote  the  usual  polar  coordinates  and  consider 


r = b(r,  0:  y) 


0 = 0 


(3.3a) 


If  we  further  require  that  there  are  no  fluctuations  in  0 , then 
the  system  in  (3.3a)  reduces  to  the  one-dimensional  equation  (3.3). 

We  assume  that  for  y > 0,  | b*  (Xg,y)  | > 0 and  that  at  y = 0,  the 

equation  exhibits  a Hopf  type  bifurcation.  Equation  (3.1)  corresponds 
to  the  stable  limit  cycle;  and  equation  (3.2)  corresponds  to  the 


l 


unstable  limit  cycle.  The  solutions  are: 


v(x)  = k exp 


b’(x0)  < 0 


u(x)  = k 


(/  i dz]  ’ 

' / exp  [-/  ft dz] dy! 


b'  (xQ)  > 0 


(3.4) 


(3.5) 


The  constants  k,  k are  determined  by  the  integrability  conditions 
and  boundary  conditions.  For  small  e,  we  use  Laplace's  method  to 
analyze  (3.4,  3.5).  We  obtain 


v(x)  ~ k exp 


“ I b ' (xQ) I (x  - xQ) 
ea (Xq) 


+ 0 (VT) 


u (x ) = k 


A 
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exp 


-b'  (xQ)  (y  - xQ) 
ea (xQ) 


dy  + 0 (v/T) 


(3.6) 


(3.7) 


Thus,  in  the  case  of  a stable  limit  cycle,  we  obtain  a locally 
Gaussian  density.  The  integral  appearing  in  (3.7)  is  the  error 
integral 


E(z) 


-1 


-s2/2, 

e ' ds  . 


(3.8) 


' -z/2 

The  correction  term  in  (3.7)  involves  E (z)  = e ' . The  error 


integral  satisfies 


d2E 


= -z 


dz 


dE 

dz 


~oo  i Zq  £ Z £ s;  oo 


(3.9) 
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The  appearance  of  Gaussian  forms  in  (3.6,  3.7)  is  due  to  the 
linearization  process  involved  in  Laplace's  method.  (Namely 


the  assumption  that  |b  (xQ)|>0  . ) 

There  are,  however,  instances  where  the  linearized  forms 
(3.6,  3.7)  break  down,  as  in  the  Hopf  bifurcation. 

When  y = 0,  i.e.,  at  the  bifurcation  value,  a one- 

dimensional Hopf  bifurcation  occurs: 


b’  (xQ,  0)  = b" ( x q , 0)  = 0 


b (xQ,  0)  t 0 


Hence,  when  applying  Laplace's  method,  for  y near  0 , we  must 
use  four  terms  in  the  Taylor  expansion  of  / b/ea  ds. 

Instead  of  the  error  integral,  we  find  (Mangel,  1977) 


x 

u(x)  = k J~  exp 


' ' 4 ' ' 3 

(b  (xQ,y)  (y-xQ)  b (xQ,y)  (y-xQ) 


6ea  (xQ) 


3ea  (Xg) 


(3.10) 


(3.11) 


% % 

where  x^x^,  x(x)  and  n(p)  are  regular  functions  of  their 
arguments  and  n(0)  = 0.  The  result  (3.12)  can  be  obtained  by 
applying  Levinson's  result  (Levinson,  1962)  directly  to  (3.5). 
Similarly,  we  find 


v(x)  'V,  c 


(3.13) 


Thus,  we  are  led  to  a new  special  function,  the  incomplete  Hopf 
integral. 


H±(z,  6) 


£ Z H 


These  integrals  satisfy 


(3.14) 


d2H+  , dH+ 

— i - ±(z  - 8z>  -Jj  (3*15) 

dz 

It  can  be  shown  that  H+(z,  6)  are  related  to  the  modified  Bessel 
functions  Kn»  (Abramowitz  and  Stegun,  1965) . It  can  also  be 

shown  that,  for  3 large,  H+(z)  ~E(z(z)),  where  z:(z)  is  a 
regular  function  of  z. 
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SECTION  4 

FIXED  LIMIT  CYCLES:  STATIONARY  ASYMPTOTIC  SOLUTIONS. 
"HINDSIGHT  AND  FORESIGHT" 

In  this  section,  we  construct  formal  asymptotic  solutions 
of  the  stationary  versions  of  (2.12,  2.13)  for  fixed  limit  cycles. 
Namely,  y is  bounded  away  from  any  bifurcation  values.  In  the 
next  section,  we  allow  y to  vary  and  consider  the  bifurcation 
case . 

UNSTABLE  LIMIT  CYCLE 

We  seek  a solution  of  the  stationary  version  of  (2.13)  in 
the  form 

u(x)  ='  £ engn(x)E(iJi  (x)//e)  + en+Jshn(x)E 
In  (4.1),  gn(x),  hn(x),  and  iMx)  are  to  be  determined.  When 

•I 

derivatives  are  evaluated,  (3.9)  is  used  to  replace  E (\p//F) 

9 

by  -E  (<J>//e)  • \p/fe~  . After  substitution  into  (4.2),  terms  are 
collected  according  to  powers  of  e.  We  obtain: 


0 = ^ e11  ^ (gn— i^hn)  ^ . tj; ) E (ij>//e) 

n=0  1 * 1 3 


n.,i  n , a1-*  n-1  . i n-l._.,  , 

+ e (b  gi  + ~2~  g.^  + c gi  )E(^//e) 


(4.2) 


' , /—  \ n+Js/.i.n  . ij  n,  . a1-1  n.  , n i, 

+ E (\ l>//e)e  hi  + a Jgi4»j  + ~y  g j + g c ipi 


i.n . . . i,n-l  ij  n,  a1-1  ,n-l  a13  .n,...  . >) 

- c h ip^  + c hi  - \pa  j~  h 


The  leading  term,  n=0  , vanishes  if 


b1ipi  - 


(4.3) 


,i  0 . 

t>  gt  = o 


(4.4) 


.i.O  . a J 0,  , ij  0,  ,0  ij  , . , 0 i,  i.O 

b hi  + — g i|>ij  + a Jg.i|;..  - h^a  + g c \p.  - c h \ 


.0 


Ji^j  - + g c ^ * c h \p\pi 


(4.5) 


-^r-hu((W.)j)  = 0 


The  transformation  <p  = -hip  converts  (4.3)  to  the  Hamilton-Jacobi 
or  eikonal  equation 


b1<|>i  + <p^<pj 


(4.6) 


I 


The  interpretation  of  the  eikonal  equation  (4.6)  in  the  stochastic 
context  is  discussed  by  Ventcel  and  Freidlin  (1970),  Ludwig (1975) , 
Mangel  and  Ludwig  (1977) , and  Mangel  (1977) . 

An  argument  using  Hamliton-Jacobi  theory  (Mangel,  1977) , 
shows  that  <p  = ip  = 0 on  the  limit  cycle  U.  Since  4>  is 
constant  on  U,  4> _ = 0 there.  We  differentiate  (4.6)  with 

U 

respect  to  x : 


b jk^i  + b^ik  + *i*j  + V ^ik^j  + Mkj*  = 


(4.7) 


Since  <p  = -Hi>  = 0 and  \p  = 0 on  U,  4>i==“ 
(4.7)  becomes 


= 0 on  U.  Thus, 


bi*ik  - 0 = §0 


(4.8) 


We  differentiate  (4.7)  with  respect  to  x1,  use  the  fact  that 
4>q  = 0 on  U and  obtain: 


ae<* 


) + 2bS_<f>  = -aSS  (<j>  ) 

ss  ,sTss  Tss 


(4.9) 


In  obtaining  (4.9),  we  have  switched  to  (s,0)  coordinates  on  U. 

(equation  (2.2)).  If  we  set  W = 4>  ^ » equation  (4.9)  becomes  a 

s s 

linear  equation  for  W: 


dW 

d9 


2bS  W 
,s 


ss 

a 


(4.10) 
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The  interpretations  of  4> , W are  interesting.  The  leading  order, 

2 

u(x)  is  constant  if  <(>(x)  is  constant.  Since  <j>(x(s))  = <J>  (6s)  /2= 

s s 

2 

(6s)  /2W,  level  curves  are  obtained  a distances  proportional  to 
1//W  . Hence  W is  a "local  variance". 

We  introduce  the  integrating  factor 


(4.11) 


and  seek  a periodic  solution  (of  period  0 ) of  (4.9).  We  obtain: 

-i  r re+©  ss  i i 

r (9)  r ass(r,e  )de 

i - (l/r  (0) ) J^e  r (9 ' ) 

Since  W(0)  is  a measure  of  the  distance  to  a contour  of  u(x), 

(4.12)  has  an  interesting  interpretation.  The  first  factor  is  a 
purely  deterministic  factor  that  causes  the  contours  to  spread 
apart.  The  second  factor  is  due  to  stochastic  effects  and  causes 
the  contours  to  close  together. 

Now  consider  (4.4).  Since  b1  = dx1/dt,  equation  (4.5) 
indicates  that  g°  is  constant  on  trajectories.  Following  Mangel 
(1977)  we  set  g°  to  be  the  same  constant  on  all  trajectories. 

This  constant  is  determined  so  that  the  leading  part  of  (4.1) 
satisfies  the  boundary  conditions.  We  set  u(x)  =0  if  x0S2 
and  u(x)  =1  if  x^S^.  Suppose  that  S^,  S2  are  level  curves 
of  ip,  with  \p  = on  S^  and  tp  = on  S2<  In  (3.9),  we  set 


(4.12) 
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Then,  to  leading  order  u=0  on  S2  and  u=l  on  S^.  If  si*S2 
are  not  level  curves  of  <|>  , we  proceed  as  follows.  Let  iJ'j'  be 
the  maximum  value  of  ip  on  and  ^ be  the  minimum  value  of 

tp  on  S2  , then  set 


(4.15) 

j 

It  can  be  shown  that,  if  ^ is  bounded  away  from  zero  on  and 

S2  then  u(x)  is  exponentially  small  on  S.^  and  l-u(x)  is 
exponentially  small  on  S2  (Mangel,  1977). 

Next,  consider  equation  (4.5).  On  U , where  <Jj=0  we  obtain 


20  = 


Z1  = 


g°  = l/EdP^/Ze) 
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Once  n is  known  on  U , it  can  be  determined  off  U by  the 
method  of  characteristics  (see  Mangel , 1977 ) . 

The  leading  part  of  the  expansion  (4.1)  is 


u(x)  ~ gUE(iJ>//e)  + 0(/e)  • 


(4.18) 


Hence,  once  g and  ip  are  known,  we  can  construct  contours  of 
u (x)  . 


STABLE  LIMIT  CYCLE 

We  now  consider  a stable  limit  cycle,  so  that  we  seek  a 
solution  of  the  stationary  version  of  (2.12): 


% (a1-^v)  . . - ( (b1  + ec1)  v)  . = 0 

■ 11  1 


Our  treatment  is  slightly  different  from  that  of  Ludwig  (1975) . 
We  seek  a Gaussian  solution  of  the  form 


v (x ) = e"^(x)  //e  ( z q + + ...  ) 


After  evaluation  of  derivatives  and  substitution  into  (4.19), 
terms  are  collected  according  to  powers  of  e (see  Ludwig,  1975) 
The  leading  term  will  vanish  if  ip  satisfies 


b^i  + j'l'  = 0 . 


(4.19) 


(4.20) 


(4.23) 


22 


The  change  in  sign  in  going  from  (4.3)  to  (4.23)  is  important 
(see  section  6)  . ir  <j>  = ^ / we  obtain  the  eikonal  equation 
(4.6),  so  that  the  analysis  of  (4.23)  is  identical  to  the  analysis 
in  the  previous  section.  We  find 


Zq  exp 

’-4>cJ6s)2  ' 

ss 

= zQ  exp 

_(6s)2 

2e 

2eW 

_ _ 

. 

(4.24) 


Hence,  eW/2  has  the  interpretation  of  a local  variance  about  the 
stable  limit  cycle.  Such  an  interpretation  has  been  given  by 
Ludwig  (1975)  . 

The  function  z^  can  be  determined  by  integration  along 
the  characteristics  of  (4.23)  (Ludwig,  1975,  Mangel,  1977).  Then 
to  leading  order 


v (x)  ^ 


zoe 


-<!'  (x)/e 


/v-*2“ 


)/e 


+ 0 ( £ ) 


(4.25) 


dx 


Ludwig  (1975)  shows  how  to  determine  z by  the  method  of  charac- 
teristics . 


SECTION  5 


HOPF  BIFURCATION 

The  analysis  of  the  preceeding  section  breaks  down  at  the 
Hopf  bifurcation,  because  the  linear  dynamics  vanish  at  the 
bifurcation  point.  The  analysis  of  section  3 suggests  a possible 
form  of  the  correct  asymptotic  solution.  The  Hopf  problem  is 
closely  related  to  a point  source  problem  for  the  wave  equation. 

Zauderer  (1970)  used  Hadamard's  method  (Hadamard,  1951)  for  such 
problems.  Our  construction  is  considerably  simpler  than  Zauderer' s, 
and  can  be  shown  to  be  equivalent  to  his. 

UNSTABLE  LIMIT  CYCLE,  STABLE  FOCUS 

We  seek  an  asymptotic  solution  of  the  stationary  backward 
equation  of  the  form 

u(x)  = ^]rngn(x)  B/e*5)  + en+i5hn(x)H  (ty/e** , B/e*5)  , (5.1) 


where  H(z,8)  = H_(z,(j)*  defined  in  (3.14).  When  derivatives  are 

" '3  3/4 

evaluated,  (3.15)  is  used  to  replace  H by  H (\p  B<^)/e 

We  assume  that  B has  an  asymptotic  expansion 


B 


(5.2) 


After  terms  are  collected  according  to  powers  of  e,  we  obtain: 
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0 = £eW  ty/e^B/e**)  ^(b1^  “ $ j 3 - 6 ) J Jgn 

n„,  . |i_i  n ^ a1"1  n_1  . i n-1 

+ e H(  T gi + — gij  +cgi  J 

+ eh+3/V  ( ) | b1hj  + cVtj  + c1i|iihh(i|,3  - eo+l 

* + ^"*ij  + "Jj1-  2hj*jV  ' 60*) 


hn  (ip  3-30<l^ ) 


(4»  - eoijo 


- hntjj_.  (tj,3  - B0<M  - hn4>  i 3 - B^j) 


n+1  fai3  / 

^ek){—  Mj(g 


n+l-k  .n+l-k  . . 3 . . 

+ h (4/  - $Qili) 


(5.3) 


+ hn+1  k(b14)i  + ip^  (ill3  - $Qip  | 


"+1  aij  r 
+ T-  "'i*jh 


- lx  *Bk(ci*ihn-k  - + hh'knj)) 

A „ ai2  . n-k ) 

]Jl6k"T^i^jh  | * 


In  (5.3),  if  a superscript  is  less  than  zero,  that  term  is  set 
equal  to  zero.  The  leading  term,  n=0  , is  composed  of  three  parts 
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First  consider  (5.4).  Since  b1  vanishes  at  the  stable  focus  P, 

3 

we  set  \p  - = 0 at  the  focus.  This  insures  that  ^ will  have 

non-vanishing  first  derivatives.  On  the  limit  cycle  U,  we  also 

3 

set  ij;  - Bq' I)  - 0.  Since  u(U)  > u(P)  we  require  that 


iMP)  =0  MU)  = /T0  • 


(5.8) 


When  the  limit  cycle  and  focus  coalesce,  we  obtain  0 = /60r 

3 

= 0.  The  sigularities  of  F (\p)  = \p  - now  match  the 
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i.e. , 


the  singularities  of  the  deterministic  system.  This  is  the  reason 
that  our  method  can  be  used  to  produce  a uniform  solution. 

The  value  of  Bg  is  still  undetermined.  It  can  be  calculated 
by  the  following  iterative  procedure.  Since  (5.4)  is  a first  order 
partial  differential  equation,  the  method  of  characteristics  can  be 
used  to  solve  it,  starting  just  off  U,  where  ^ = /s  ^ and  8 ^ is  the 
initial  estimate  for  0q.  We  follow  characteristics  (.called 
rays)  that  approach  P.  If  ip  does  not  approach  0,  then  0 ^ 
must  be  replaced  by  a better  estimate  B^1^  . The  method  of  false 

position  can  be  used  to  calculate  iterates  of  Bg*  In  this  fashion, 

Bg  can  be  determined  to  any  order  of  accuracy.  An  alternative 
procedure  would  follow  rays  from  P to  U.  The  choice  of  method 
must  be  made  on  the  basis  of  numerical  practicality. 

Although  (5.4)  can  be  solved  by  the  method  of  characteristics, 
our  main  interest  is  in  experiments  beginning  near  U.  Consequently, 
we  determine  ip  in  a vicinity  of  the  limit  cycle  by  a Taylor  ex- 
pansion. We  assume  that  B > 0.  Equation  (5.4)  is  differentiated 
with  respect  to  x and  then  changed  to  (s,0)  coordinates.  We 
obtain: 

dij;  , 

-gf  + b®g(0,  V>H»S  - aVg  = 0 • (5.9) 

In  deriving  (5.9),  we  have  used  the  fact  that  ^ = /8Q  on  U (so 
that  3ip2  - Bg  = 2Bq  on  U)  . 
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Equation  (5.9)  is  a version  of  Abel's  equation  (Davis,  1962).  We 
introduce  a new  variable  z,  defined  by 


1 p = 1/Bz 


where  B = b B 
*s 


The  periodic  solution  of  (5.9)  is  then 


4'  (0,  9,  y ) = 


O-t-fc) 

-2  / 51 

•'9  B( 


bX  I4fr  ds 


s>  ds  . 

s)  r (©)  - l 


where 


(5.10) 


(5.11) 


B ( s ) = exp 


b*  (0,  9,  y ) d0 

9 S 


(5.12) 


r (s)  = 1/b  (s) 


(5.13) 


Equation  (5.5)  indicates  the  g®  is  a constant.  The  value  of 
can  be  determined  exactly  as  in  section  4.  Equation  (5.6)  is 
analogous  to  (4.5).  It  is  slightly  more  complicated  since  it 
contains  the  unknown  paramter  6-^.  This  parameter  can  be  determined 
in  the  same  fashion  as  was  determined. 

It  can  be  shown  that  all  of  these  constrictions  are  regular 
at  the  bifurcation  point  y = 0 . The  proof  is  analygous  to  the 


proof  given  in  Mangel  (1977)  for  other  stochastic  dynamical  systems. 
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STABLE  LIMIT  CYCLE,  UNSTABLE  FOCUS 


In  this  case,  we  are  interested  in  uniform  solutions  of  the 
forward  equation  (4.9).  It  is  clear  that  the  Gaussin  ansatz  in 
section  4.2  breaks  down  for  y small.  We  seek  a solution  of  the 


form 


v(x)  ~ exp 


f I / <J>  (x  )4  (x  )2\ ' 


(Z°(X  ) + E Z 1 (x ) + ...)  I 


Following  the  procedure  of  section  4.2,  we  are  led  to 


14'i  + — 2~  ♦ - 80^)  = 0 


The  change  is  sine  in  going  from  (5.4)  to  (5.19)  is  important. 
Equation  (5.19)  can  be  treated  by  the  method  of  characteristics  or 
by  a Taylor  expansion.  The  function  z^(x)  can  be  determined  by 
integration  along  the  characteristics  of  5.15  (Ludwig,  1975). 

Thus,  the  stationary  distributions  for  the  Hopf  bifurcation 
problem  have  been  determined.  These  distributions  are  regular 
functions  of  vi,  the  deterministic  bifurcation  parameter. 


5.14) 


(5.15) 


SECTION  6 


FIXED  LIMIT  CYCLES:  SOME  EXAMPLES 

In  this  section,  we  present  a number  of  numerical  examples 
that  illustrate  the  behavior  of  u(x)  and  v(x),  as  determined 
in  section  4.  For  convenience,  we  use  systems  already  in  polar 
coordinates. 


(6.1) 
(6.2) 

(6.3) 

The  circle  r = 1 is  an  unstable  limit  cycle.  Let 

u(x)  = Pr{process  hits  r = 1.98  before  r = . 02|x(0)  = x} 

In  figure  3,  we  show  the  u = .8,  .9  contours  6r(0)  where  0 
measures  distance  along  the  cycle  and  6r  is  the  distance  from 
r = 1 to  the  contour.  The  noise  and  deteministic  dynamics  are  in 

I 

phase.  Both  contours  are  of  the  form  6r(0)  = k (k  + cos0)  where 

I » 

k,  k are  constants  and  k >1.  in  table  1,  we  compare  the  theory 
with  Monte  Carlo  experiments. 


EXAMPLE  6 . 1 

r = r (r-1)  (2-r)  (1 . 1 + cos0) 

0 = 1 

with  covariance 

earr  = .l[r2  + (2-r)?J(1.5  + cos0)2 


I 


TABLE  1 

COMPARISON  OF  THEORETICAL  AND  MONTE  CARLO  RESULTS 

FOR  EXAMPLE  6.1 


Initial  point 

u (Theory) 

u (Monte  Carlo)* 

(1.10, 

0) 

.60 

.57 

(1.21, 

.21) 

.70 

.68 

(1.50, 

.42) 

.90 

.92 

(1.30, 

3.97) 

.80 

.79 

(1.53, 

5.86) 

.90 

.92 

(1.08, 

1.88) 

.60 

.58 

(1.19, 

2.30) 

.70 

.69 
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EXAMPLE  6.2 

We  now  consider  the  system 

r = r (r-1)  (r-2)  (1. 1 + cos6)  (6.5) 

6=1  (6.6) 

with  earr  given  by  (6.3).  The  deterministic  system  has  a stable 
limit  cycle  at  r = 1.  Let 

v(x)dx  = Pr{process  is  eventually  found  between  (x,  x + dx) } (6.7) 

In  figure  4,  we  plot  the  .91,  .99  contours  of  v(x)  as  a function  of 
6r(0)  , where  6r  is  the  distance  from  the  cycle  to  the  contour. 

EXAMPLE  6.3 

We  now  take 

r = r (r-1)  (2-r)  (1. 1 + sin0)  (6.8) 

rr 

with  ea  given  by  (6.3).  In  this  case,  the  noise  is  out  of 
phase  with  the  deterministic  cycling.  In  figure  5,  we  plot  the 
u = .8,  .9  contours  and  in  table  2,  compare  Monte  Carlo  and  theore- 
tical results. 

If  we  take 

r = r (r-1) (r-2) (1. 1 + sin0)  (6.9) 
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TABLE  2 


COMPARISON  OF  THEORETICAL  AND  MONTE  CARLO  RESULTS 

FOR  EXAMPLE  6.3 


Initial  Point 

u (Theory) 

u (Monte  Carlo)* 

(1.14,  0) 

.60 

.62 

(1.22,  .21) 

.70 

.73 

(1.49,  .42) 

.90 

.90 

(1.55,  3.97) 

.80 

.82 

(1.66,  5.86) 

.90 

.87 

(1.12,  1.88) 

.70 

.67 

(1.06,  1.47) 

.60 

.64 

then  r = 1 is  a stable  limit  cycle.  In  figure  6,  we  plot  the 

v = .91,  .99  contours.  We  note  the  shift  in  phase  of  u,v. 

to  r 

For  fixed  r , ea  reaches  its  maximum  at  0,211.  The  behavior 

of  the  solution  of  the  backward  equation,  u(x),  "anticipates"  the 

• rr 

noise,  in  that  6r(0)  reaches  its  maximum  before  ea  (0)  reaches 
its  maximum.  The  density  v(x)  solution  of  the  forward  equation, 
exhibits  hindsight  in  the  6r(0)  peaks  after  the  maximum  value  of 
the  noise. 

Comparison  of  figures  3 and  5 is  also  interesting,  in  light 
of  the  interpretation  given  to  W in  section  4.  Namely,  the 
deterministic  and  stochastic  terms  in  (4.12)  "compete"  with  each 
other,  the  former  increasing  W and  the  latter  decreasing  W. 

For  the  situation  exhibited  in  figure  3,  the  deterministic  and 
stochastic  terms  are  "in  phase"  and  the  contours  are  relatively 
constant.  On  the  other  had,  for  the  situation  exhibited  in  figure 
5,  the  deterministic  and  stochastic  terms  are  out  of  phase.  The 
contours  exhibit  sinusoidal  oscillations. 
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